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The fully nonlinear notion of resonance— geometrical resonance— in the general context of 
dissipative systems subjected to nonsteady potentials is discussed. It is demonstrated that 
there is an exact local invariant associated with each geometrical resonance solution which 
reduces to the system’s energy when the potential is steady. The geometrical resonance so¬ 
lutions represent a local symmetry whose critical breaking leads to a new analytical criterion 
for the order-chaos threshold. This physical criterion is deduced in the co-moving frame from 
the local energy conservation over the shortest significant timescale. Remarkably, the new 
criterion for the onset of chaos is shown to be valid over large regions of parameter space, 
thus being useful beyond the perturbative regime and the scope of current mathematical 
techniques. 

Hamiltonian and dissipative systems have traditionally been studied separately due to their 
clearly different dynamic properties 1 : dissipation forces give rise to the existence of transient dy¬ 
namics associated with the basins of the different attractors, while the Poincare integral invariants 
of Hamiltonian systems lead to special behaviour of the eigenvalues of equilibria and periodic orbits, 
and to existence theorems for various types of orbits such as the celebrated Kolmogorov-Arnold- 
Moser theorem. To date, only the notion of geometrical resonance (GR) 2 has been able to provide 
a deep link between autonomous Hamiltonian and non-autonomous dissipative systems in the sense 
that it offers a universal procedure with which to locally ” Hamiltonianize” an otherwise dissipative 
system by suitably choosing the non-autonomous term(s) /*(t) such that the system’s energy is con¬ 
served locally: /*(t) = The original formulation of GR analysis was for steady potentials 2 , 

and was applied to diverse nonlinear problems involving such potentials 3-11 , including the stability 
of the responses of an overdamped bistable system', the suppression of spatio-temporal chaos and 
the stabilization of localized solutions in general spatio-temporal systems 5,8-10 , quantum control 
in trapped Bose-Einstein condensates (BECs) 11 , and the characterization of periodic solutions of 
a fractional Duffing’s equation 12 . 
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On the other hand, a major body of research has considered nonsteady potentials appearing 
in different physical contexts such as synchrotron motion of beams 12 , BECs in optical lattices 14 , 
electron transport in semiconductor superlattices 15 , and nanoscale devices powered by the lat¬ 
eral Casimir force 16 , just to cite a few representative examples. In general, the reference frames 
co-moving with nonsteady potentials are accelerating frames, which introduces an additional com¬ 
plexity into analysis of the dynamics relative to the laboratory reference frame (L-frame). Since 
GR is neither more nor less than a local symmetry, namely that the dynamics equations remain 
locally invariant under time reversal when the non-autonomous terms are suitably locally chosen, 
fi(t) = fi.cmit), it seems appropriate and pertinent to explore its implications in general systems. 

In this paper, new properties of this subtle symmetry in the generalized context of dissipative 
systems in nonsteady potentials are characterized and exploited to determine a physical criterion 
for the order-chaos threshold in parameter space whose accuracy and scope goes beyond current 
perturbative mathematical techniques. The theory is discussed through the universal model of 
a particle subjected to a spatially periodic and temporally shaken potential which describes, for 
example, the chaotic phase oscillation of a proton beam in a cooler synchrotron 13 . 

Results 

Theory 

Let us consider the class of non-autonomous dissipative systems 

x = g{x,t) - d(x,x) , (1) 

where the overdot denotes d/dt, g(x,t ) = — dV (x,t) /dx, with V (x,t) = V [x — f(t )] being a 
nonsteady and spatially periodic potential while f (t) is an a priori arbitrary function of time, 
and where — d (x, x) is a generic dissipative force. In the potential reference frame (V-frame) with 
z (t) = x ( t) — f(t), Eq. (1) reads 

z = 9{z)~ d(z,z,fj,tj - f. (2) 

In general, if zanit) is a GR solution of Eq. (2), it must satisfy 

ZGR = g ( ZGR ) , (3) 

fGR=- d ( z GR,ZGR,fGR,fGR, t )i ( 4 ) 

and hence 


XGR = ZGR. + fcR 


( 5 ) 
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is a GR solution of Eq. (1) for a given set of initial conditions [x'o = x(t = 0), xq = x (t = 0)]. 
Definitions (3)-(5) give rise to the following distinguishing properties. First, in contrast with 
the case of steady potentials 2 , where f GR is univocally determined from an algebraic equation 
involving the single GR solution associated with a given set of initial conditions, one has to solve a 
differential equation for f G R in the present general case and obtain the initial values fen (t = 0) = 
Igr, o, fcR(t = 0) = / gr,o as a P ar t °f the whole solution. This is because the GR scenario for 
nonsteady potentials involves two reference frames, the V-frame being non-inertial in the general 
case. Second, conditions (3)-(4) are equivalent to the local energy conservation requirement 

(1/2 )z GR + V [ z G r} = E gr = const (6) 

in the V-frame, while one has the requirement of a different local invariant in the L-frame: 

- xcR^j = const. (7) 

After Taylor expanding the potential, the local invariant I GR can be recast into the more trans¬ 
parent form 

(l/2)x GR + V [x GR - f G R,o] + S [x GR , f GR } = const , (8) 


Igr = 


_ X GR 


+ V [ x gr ~ Igr] + Igr 


where, 

S [ x gr, Igr] = ^ ~J~(fGR ~ fGR,o) n V (n) [x GR - /gr,o] + f G r (/ g ^/ 2 - Igr) , ( 9 ) 

71=1 

with V ( n ) = d n V(z)/dz n . From Eq. (8) one sees that, under GR conditions, the energy associated 
with the corresponding steady potential is not (locally) conserved in the L-frame, as expected, 
while the new invariant I G r allows the temporal evolution of this energy to be calculated for each 
GR solution. Third, in the Hamiltonian limiting case, i.e., d (x, i:) —>• 0, Eq. (4) becomes f GR = 0, 
and hence f G n{t ) = Vt with V being an arbitrary constant and where an additional additive 
constant has been taken to be zero without loss of generality. This means that, in the absence 
of dissipation, GR solutions are solely possible for potentials traveling with constant speed, i.e., 
for inertial frames, as expected. And fourth, a GR solution will be observed if it is stable, i.e., if 
any small perturbation 8z ( t ) of z G r ( t) is damped. After substituting z(t) = z G r ( t ) + 8z ( t) into 
Eq. (2) with / ( t ) = f G R ( t ), one obtains the linearized equation of motion for small perturbations 
5z ( t ): 


Sz + 



Sz + 
GR 





( 10 ) 
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Note that this generalized Hill equation with dissipation also governs the stability of the GR 
solutions in the L-frame (cf. Eq. (5)). It is shown below that this stability analysis together with 
the dependence of the GR solutions on the system’s parameters and the local invariants (6) and (8) 
allows one to get a new analytical criterion for the order-chaos threshold from the weakest useful 
approximation to the local energy conservation in the V-frame. 

Universal model 

To demonstrate the effectiveness of the present GR theory in a simple paradigmatic model, 
consider the dissipative dynamics of a particle subjected to a spatially periodic and temporally 
shaken potential: 


x + sin [x — f(t)] = — rjx. (11) 

The dimensionless Eq. (11) describes for example the pinion motion of a nanoscale device composed 
of a pinion and a rack coupled via the lateral Casimir force, where i] > 0 is a damping coefficient 
while f(t) accounts for the a priori arbitrary horizontal motion of the rack 16 . In the V-frame with 
z (t) = x (t) — f(t), Eq. (11) reads 

z + sin z = —rjz — rjf — f. (12) 

Thus, GR solutions of Eq. (12) must satisfy 

z G r + sin z G r = 0, (13) 

zgr + fcR + Igr/v = 0- (14) 

Exact analytical periodic solutions of the integrable pendulum (13) 17 corresponding to libration 
and rotation motions are given by 

z l GR (t]t 0 ,m) = 2arcsin [y/msn (t — to; m)l , 

z G r( t', to,m) = 2\/mcn (t — to] m ), (15) 

and 


ZG R {t] to, m) = ±2 am [(t — to) /\fm.] m , 


z r G R {t]to,m) = ±—)= dn 


m 


t -1 0 


; m 


m 


(16) 


respectively, where sn(-;m), cn(-;m), dn(-;m), am(-;m) are Jacobian elliptic functions of pa¬ 
rameter m € ]0,1[, the upper (lower) sign in the rotation solutions refers to counterclockwise 
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(clockwise) rotations, while to is an arbitrary initial time. These solutions have the respective 
periods T l (m) = 4 K and T r (m ) = 2y/mK, where K = K(m ) is the complete elliptic integral 
of the first kind 18 . Although the parameters corresponding to libration and rotation motions are 
inversely related each other, the same notation, m, is used here since both parameters are defined 
over the same interval 17 and there is no possibility of confusion in the subsequent analysis. Also, 
definition sn(-;m) = sin[am(-;m)] has been used to write a simpler alternative expression for 
z Q R (t\ tQ.m). After taking to = 0 for simplicity, using the Fourier series of the Jacobian elliptic 
functions involved 18 , and integrating the linear differential equation (14), one straightforwardly 
obtains the GR excitations 


Igr( i) = C*i + C 2 e b n cos ( u> n t + ip n ), 

?i=0 

TTt °° 

fcR^f) = c[ + C' 2 e r,t =F ^ X] ^ n cos + Pn) > (1?) 

* n= 1 

and the corresponding GR solutions in the L-frame 

x b 

x Gr( t) = Cl + C *2 e r,t + 7-sin ( w nt + Pn) > 

n= 0^ ^ n 
o° , 

*&*(*) = C[ + C' 2 e~ vt shl fat + vQ , (18) 

n= 1 S ^ n 

with 


bn — 
6( = 


27rr/ 


: sech 


UnK^Jr] 2 + uj 2 
2irri 

u' n Ky/m y 7 i/ 2 + t<4 2 
= (n+ 1/2) 7 t/A', 

w n = W [VrnK) , 

Pn = arctan (rj/u n ), 

^4 = arctan ( tj / lj ' n ) 


(2n + 1) 7 t/G 


sech 


2K 

( mrK' 


V 


K 


(19) 


and where C'i ; 2 ,C'( 2 are constants to be determined from the initial conditions (xo,io) (see the 
Methods section), K' = it (1 — m), while the upper (lower) sign in Eqs. (17) and (18) refers to 
counterclockwise (clockwise) rotations. These GR solutions have the following properties, (i) Their 
stability is governed by Eq. (10), i.e., 

Sz + rj5z + cos Z GR$ Z ~ 


( 20 ) 



6 


which reduces to the Lame equations 

d 2 u 

+ [l — f? 2 /4 — 2m sn 2 (t; m)] u = 0, (21) 

d 2/ u 

-pj + \m (l — ? 7 2 / 4 ) — 2m sn 2 (r; m)] v = 0, (22) 

where u = exp (rjt/2) 8z and v = exp (jii/mr/ 2) 5z for librations and rotations, respectively. Stan¬ 
dard results for these Lame equations 19,20 indicate that Eq. (20) presents only one instability region 
in the m — if parameter plane. A careful comparison of Eq. (21) with Eq. (22) leads one to ex¬ 
pect the instability region for librations to be clearly narrower than that for rotations owing to 
the term m (1 — 7/ 2 /4) < 1 — rj 2 /4 since m € ]0,1[. Moreover, the maximum range of rj values 
in the instability regions is expected to occur when m ~ 1 for both kinds of motion due to all 
GR solutions z 1 q R (t) converging to the separatrix (the most unstable phase path) of the integrable 
pendulum asm-} 1. Numerical simulations confirmed these expectations, as is shown in Fig. 1. 
(ii) For any set of initial conditions not on the unperturbed separatrix, i.e., for any GR excita¬ 
tion (17) and corresponding solution (18), one sees that the dependence of each harmonic of such 
excitations and solutions on the damping coefficient has exactly the same form: rj/y/i] 2 + a, with 
a being a function of the corresponding natural period. From this it can be inferred that, for a 
periodic excitation f(t) of amplitude 7 , the dependence of the chaotic-threshold amplitude, 7 th, on 
r/ should obey this functional form irrespective of the value of 7 , an unanticipated result in view 
of the perturbative character of the current mathematical techniques to predict the onset of chaos 
(Melnikov’s method (MM)), (iii) Since harmonic functions are commonly used to model periodic 
excitations, the GR solution corresponding to libration near the bottom of the potential well (i.e., 
xo,~ 0,io ~ 0, m > 0) is of especial interest. One straightforwardly obtains the steady (t 3 > 7]~ 3 ) 
solutions (cf. Eqs. (17) and (18)) 

x l G R(t) = 2 y / msini + f l GR {t), 

f l GR (t) = cos [t -g arctan 7 ) + O , (23) 

V1 + V 2 V 2 

while the corresponding local invariant ( 8 ) reduces to 

ll GR = 0 -/ 2 ) x gr (t) - cos X GR (t) + - GR — * + O (m 3 ) , (24) 

i.e., I l GR (t) is no more than the sum of the energy associated with the limiting case of steady 
potential plus the energy associated with the V-frame moving as a linear harmonic oscillator of 
period 2ir , which is an unexpected result, (iv) In the Hamiltonian limiting case, i.e., i] —>• 0, GR 
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solutions for librations are not possible due to their oscillatory character around a fixed point 
is incompatible with the requirement of a traveling potential function (cf. third property in the 
previous subsection), while GR solutions for rotations are indeed possible (see Eq. (17)). (v) In 
the limit of very high dissipation (77 —>• 00 ), the steady GR solutions are equilibria (cf. Eqs. (18) 
and (23)), as expected. 



FIG. 1: Stability-instability charts. Stability-instability charts obtained by numerical integration of 
the Lame’s equations (a) Eq. (21) and (b) Eq. (22) for libration and rotation motions, respectively, where 
instability regions are indicated by dots. The insets show the tips of the instability tongues when m —> 1 
and 77 —>• 0 . 

Order-chaos threshold 

Next, one can use the above properties of the GR solutions to obtain an analytical estimate of 
the order-chaos threshold associated with a generic T-periodic excitation f g (t) of amplitude 7 . In 
general this generic excitation will not exactly correspond to any GR excitation function (17), and 
hence one cannot expect strict conservation (i.e., over an infinite timescale) of the invariants ( 6 ) 
and (7) for any set of initial conditions. Indeed, the energy flow is governed in the V-frame by the 
equation (cf. Eq. (12)) 

^ = —rjz 2 - nzf g - zf g . (25) 

For each set of initial conditions, the closer the excitation f g (t) is to the corresponding fcR(t), 
the smaller the deviation of the energy from the corresponding local invariant Eqr- Clearly, the 
weakest physical condition that will cope with this deviation is that the energy be locally conserved 
over the shortest significant timescale, i.e., as an average over a period of the corresponding GR 






















solution: 


(^jf) Tl ^ = j Tl r (p Z GR + ll z GRfg + z GRfg) dt = 0 ( 26 ) 

for some to, and where z 1 q R = z 1 q R (f;to,m). Also, one assumes a Galilean resonance condition—a 
necessary condition for GR (cf. Eq. (14))—for both libration (T = T l /(2n + l) for some n = 0,1,...) 
and rotation (T = T r /n for some n = 1,2,...) motions. Thus, Eq. (26) provides a local condition 
that takes into account the initial phase difference between the generic excitation and the GR 
solution, hence allowing one to obtain a threshold condition (in particular, a threshold amplitude 
- jth ) for the energy conservation in its weakest sense. According to the above stability analysis, GR 
solutions are not uniformly stable as the natural period is varied. Therefore, Eq. (26) is subject 
to the caveat that it is not expected to be uniformly valid for all values of the excitation period 
because of its dependence on the integration domain. In the limiting case T l,r —>• oo, when both 
libration and rotation GR solutions converge to the separatrix 

z s ,±(t] to) = ± arctan [sinh [t — to)], 

z 3 ,± (t',t 0 ) = ±2sech(t - t 0 ), (27) 

the corresponding GR excitation is no longer a periodic function, as expected, but (for t 0 = 0) 

nihil 1 ± 2 -3±!Z-_e 2t j 

fo R (t ) = Cr + C 2 e~ vt =f 2 gd (t) ± 4 e*- ’ ^ ^ ’ --, (28) 

where C\p ar e constants to be determined from the initial conditions (.To,xo) (see the Methods 
section), gd (t) and 2 -Ei(E — e 2t ) are the Gudermannian and the hypergeometric functions, 

respectively 22 , and Eq. (26) becomes 

\~dt) J yl z2 s,± + ‘'l z s,±fg + z s,±fg)dt = ° (29) 

for some to- Since the separatrix is the most unstable phase path, in the sense that it is the 
boundary between two distinctly different types of motions, one would expect the onset of chaos 
when a gradual breaking of the GR local symmetry reaches a critical value. Indeed, Eq. (29) 
provides the physical condition for such a critical breaking, hence allowing the order-chaos threshold 
in parameter space to be estimated. It should be stressed that, because the GR local symmetry 
is defined over the complete parameter space, condition (29) is postulated irrespective of the 
parameter values. For the sake of clarity, consider the application of condition (26) to the simple 
case of a harmonic excitation /(f) = 7 cos (2nt/T). After some simple algebra (see the Methods 
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section), one straightforwardly obtains the following threshold amplitudes from Eq. (26): 

4 rj [E — (1 — m)K\ 


1th = 


1th = 


•Kojy/rf + cj 2 sech ( uK') 
Arj^/rriE 


(30) 

(31) 


rf + uj 2 sech ( Uy/mK') 

where oj = 2tt/T and E = E(m) is the complete elliptic integral of the second kind 18 . Also, 
limm^i 1th = 1th with 

4 rf 


'Ith 


(32) 


ttu y 7 r / 2 + oj 2 sech (ttuj/2) 
being the explicit estimate of the order-chaos threshold in parameter space. 

Let us now compare the prediction (32) with that obtained from MM. In keeping with the 
assumptions of the MM 23,1 ’ 21 , here it is assumed that one can write rj = efj, 7 = £7 where 0 < e <C 1 
and 77 , 7 , 0 ; are of order unity. Next, one calculates the Melnikov function (MF), M(to), for the 
system ( 12 ) with the harmonic excitation: 


z + sin 2 : = —erjz + eycn 2 cos (u;t) + O (e 2 ) 


(33) 


Since the MF provides an O (e) estimate of the distance between the stable and unstable manifolds 
of the perturbed system in the Poincare section at to, one readily obtains 


M{to) = —8?? ± 27T70 ; 2 sech (7ru;/2) cos (wfo) • 


(34) 


If the MF has a simple zero, then a homoclinic bifurcation occurs, signifying the appearance of 
homoclinic chaos 24 . This yields the threshold value 


ryM . 

1th 


477 


7tu ; 2 sech (irbo/2) 


(35) 


Numerical simulations confirmed the effectiveness of estimate (32) as against (35). An illustrative 
example is shown in Fig. 2, in which one sees how the chaotic regions in the 7 — 7 parameter plane, 
determined by Lyapunov exponent (LE) calculations (see the Methods section), are reasonably 
well bounded by estimate (32), while the extrapolation (recall that 77,7 must be much smaller 
than unity) of the MM estimate (35) clearly fails. Note that estimates (32) and (35) coincide for 
sufficiently small values of 77 (perturbative regime). This is not so surprising since the MF is, up to 
a constant, exactly the integral that Poincare derived from Hamilton-Jacobi theory to obtain his 
celebrated obstruction to integrability 27 , while the GR local symmetry implies a local restoration 
of integrability of an otherwise non-integrable system. Figure 3 shows illustrative examples of the 
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regularization routes as 7 and ij are changed while crossing the order-chaos threshold. Typically, 
the system ( 11 ) goes from a period -1 attractor to a strange chaotic attractor as the excitation 
amplitude increases for a sufficiently small value of the damping coefficient (see Fig. 3(a)). The 
overall evolution of the initial periodic state is characterized by the energy E = x /2 + 1 — cos x 
undergoing a period-doubling route as 7 is increased. Also, for fixed 7, the system ( 11 ) goes from 
the strange chaotic attractor existing at a sufficiently small value of the damping coefficient to a 
period-1 attractor as 77 is increased via an inverse period-doubling route (see Fig. 3(b)). 

Discussion 

In conclusion, a theory of geometrical resonance in dissipative systems subjected to nonsteady 
potentials has been presented, and its effectiveness in obtaining an analytical criterion for the order- 
chaos threshold in parameter space beyond the perturbative regime demonstrated by means of a 
paradigmatic example. From a theoretical point of view, the characterization and determination of 
the frontiers between chaotic and regular motions of real-world systems is a general problem that 
needs to be addressed in all branches of science. While the mathematical theory of deterministic 
chaos was definitively established by the work of Poincare, Birkhoff, and Smale, the present physical 
theory suggests understanding the onset of homoclinic chaos as being coincident with a critical 
breaking of the geometrical resonance local symmetry , specifically, as coinciding with a critical 
breaking of the local conservation of the separatrix’s energy in the co-moving frame for generic 
(periodic) nonsteady potentials. Finally, a natural continuation of the present work is the study 
of the geometrical resonance local symmetry and its eventual breakage in classical and quantum 
Hamiltonian systems. 

Methods 

Determination of the initial value problem. From Eqs. (13), (15), (16) with to = 0, one 
obtains 


t'cRA) — fGR h — 0) — x GR h — 0) Z GR ^ — 0) — Xo, 

f GR,o = fGR& = 0) = x GR (t = 0) — z GR (t = 0) = xo — 2 y/m, 

f GR,0 = f GR^f = 0) = X GR (t = 0) — Z GRif = 0) = x 0 =F ~/=, 

y/m 


while from Eq. (18) one obtains 


Ci 


C 2 


xo + 


Xo 


E 


b n U n COS 2 ip n 


sm (f n 

b n u n COS 2 (fin 


n =0 

OO k , , ™ 0 2 

*o + £ 

n =0 


OO 

- COS (fn, 

n =0 


(36) 


Sill (f n 


(37) 
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(b) 




(c) 



FIG. 2: Chaotic regions in the control parameter plane. Chaotic regions (dots) in the r] — 7 parameter 
plane corresponding to Eq. (11) with /(f) = 7 cos (2nt/T) for three values of the driving period: (a) T = 7.5, 
(b) T = 9.38, (c) T = 10.5. A dot is plotted on a 100 x 140 grid when the corresponding maximal LE is 
greater than 10 -3 . Dashed and solid lines represent the theoretical chaotic thresholds [cf. Eqs.(32) and (35), 
respectively] from GR analysis and MM, respectively. 


for libration motions, and 


C'i — Xq + 


1 




b' n u' n COS 2 (p' n 


C ' 2 = - 
V 


V 

-x 0 ±J2 

71=1 


sin ip' n 

b' n u'n COS 2 <p' n 


n= 1 

00 7/ , 2 




n HO® Tbi; 


n= 1 


sm 


(38) 


for rotation motions, and where b n , u n , p n , b' n ,uj ' n , <p' n are given by Eq. (19), while the upper (lower) 
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FIG. 3: Bifurcation diagrams of Eq. (11) with f(t) = 700 s (27rt/T). (a) Bifurcation diagram of 
energy E = x /2 + 1 — cos a; as a function of the amplitude 7 for 77 = 0.95. The vertical dashed line 
indicate the threshold amplitude 7 th = 2.4957 predicted from Eq. (32). (b) Bifurcation diagram of energy 
E = x /2 + 1 — cos 2 as a function of the damping coefficient 77 for 7 = 2.4. The vertical dashed line indicate 
the threshold value of the damping coefficient i~i t h = 0.8514 predicted from Eq. (32). Fixed parameter: 
T = 9.38. 


sign in Eqs. (36) and (38) refers to counterclockwise (clockwise) rotations. 

Derivation of the geometrical resonance excitation for the separatrix. For the sepa- 
ratrix (27), Eq. (14) reduces to the linear differential equation 

Ten + vIgr = =F2 rj sech (t ). (39) 


After using the method of variation of parameters 25 , one straightforwardly obtains the general 
solution given by Eq. (28) while its derivative is written 


fGR^) = -C 2 Ve vt = F 4/7 e' 


1+77. 3+77. _ 2 1 


t 2^1(1, 2 > 2 ! 


-e 2t ) 


I + 77 


(40) 


Finally, after taking into account (36) for m = 1, i.e., fQ R (t = 0) = xq, f GR (t = 0) = xq =F 2, the 
integration constants are given by, 




C 2 


4— (xo =F 2) , 
rj 


1 

V 


-x 0 ± 2 =F V 






(41) 


where if) (rj) is the psi (Digamma) function 22 , while the upper (lower) sign in Eqs. (28), (39), (40), 
and (41) refers to counterclockwise (clockwise) rotations. 

Derivation of the threshold amplitudes for the onset of chaos. Let us consider the 
simple case of a harmonic excitation f(t) = ycos (2nt/T). 
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Libration motions. Equation (26) for T l ( m ) = 4 K and the Galilean resonance condition T = 
T l / (2 n + 1), n = 0,1,reduces to 
'(1E\ r2K 


dt / 


T l 


(to) - 


' —2K 


dz l GR ~ Tt uz G r sin M) - 7 u 2 z gr cos (wt) 


dt. 


(42) 


After substituting z GR (t;to,m ) from Eq. (15) into Eq. (42), using the Fourier series of the elliptic 
function cn(-;m) 18 , and using standard tables of integrals 28 , one obtains the average energy over 
the period T l as a function of the initial phase difference (coto) between the harmonic excitation 
and the GR solution 


> (to) = 16?? [E — (1 — m)K] — sech [coK'] sin (uito) — An'-fu 2 sech [wE'] cos (coto) • 

IT 1 

(43) 

From Eq. (43) one sees that a necessary condition for {^) T i (to) to change sign at some to is 
written 



I 677 [E — (1 — m)K] ^ yf(Anjurj sech [ujK 1 ]) 2 + {Att'juj 2 sech [co/G]) 2 , 


(44) 


where the equals sign in Eq. (44) yields the threshold amplitude 7 ^ given by Eq. (30). 

Rotation motions. Equation (26) for T r ( m ) = 2\fmK and the Galilean resonance condition 
T = T r /n, n = 1,2,..., reduces to 


dE\ 

dt j rp r 


(to) - 


ry/rnK f 
' —y/mK 


VZgr ~ rfwz r GR sin (cut) - 7 u z z GR cos (cot) 


2 • r 


dt. 


(45) 


After substituting ZQ R (t-,to,m) from Eq. (16) into Eq. (45), using the Fourier series of the elliptic 
function dn(-;m) 18 , and using standard tables of integrals 28 , one obtains the average energy over 
the period T r as a function of the initial phase difference (coto) between the harmonic excitation 
and the GR solution 


(to) = 8rjy/rnE =F 2 ^ 7770 ; sech [uy/mK'] sin (coto) =F 2'k^/uj 2 sech [c oy/mK'] cos (coto) j (46) 

'j'r 

where the upper (lower) sign in Eq. (46) refers to counterclockwise (clockwise) rotations. From 
Eq. (46) one sees that a necessary condition for (y|^) Tr (to) to change sign at some to is written 

8t]^/rnE ^ \J ( 2 ^ 7 co?/sech c osfmK' ) 2 + (27T7co 2 sech ujy/mK' ) 2 , (47) 



where the equals sign in Eq. (47) yields the threshold amplitude j[ h given by Eq. (31). 

Numerical simulations. Equations (11), (21), and (22) have been numerically solved using 
a fourth Runge-Kutta method with discrete time step St = 0.001. Lyapunov exponents have been 
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computed using a version of the algorithm introduced in reference 26, with integration typically 
up to 2 x 10 4 drive cycles for each fixed set of parameters. 
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